home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / sfit.pro < prev    next >
Text File  |  1997-07-08  |  2KB  |  67 lines
  1. ; $Id: sfit.pro,v 1.4 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1993-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5.  
  6. function sfit, z, degree, kx=kx
  7. ;+
  8. ; NAME:
  9. ;    SFIT
  10. ;
  11. ; PURPOSE:
  12. ;    This function determines a polynomial fit to a surface.
  13. ;
  14. ; CATEGORY:
  15. ;    Curve and surface fitting.
  16. ;
  17. ; CALLING SEQUENCE:
  18. ;    Result = SFIT(Data, Degree)
  19. ;
  20. ; INPUTS:
  21. ;     Data:    The two-dimensional array of data to fit. The sizes of
  22. ;        the dimensions may be unequal.
  23. ;
  24. ;    Degree:    The maximum degree of fit (in one dimension).
  25. ;
  26. ; OUTPUT:
  27. ;    This function returns a fitted array.
  28. ;
  29. ; OUTPUT KEYWORDS:
  30. ;    Kx:    The array of coefficients for a polynomial function
  31. ;        of x and y to fit data.
  32. ;        This parameter is returned as a (Degree+1) by (Degree+1) 
  33. ;        element array.
  34. ;
  35. ; PROCEDURE:
  36. ;     Fit a 2D array Z as a polynomial function of x and y.
  37. ;     The function fitted is:
  38. ;          F(x,y) = Sum over i and j of kx(j,i) * x^i * y^j
  39. ;     where kx is returned as a keyword.
  40. ;
  41. ; MODIFICATION HISTORY:
  42. ;    July, 1993, DMS        Initial creation
  43. ;
  44. ;-
  45.  
  46.    on_error, 2
  47.  
  48.    s = size(z)
  49.    nx = s[1]
  50.    ny = s[2]
  51.    m = nx * ny        ;# of points to fit
  52.    n2=(degree+1)^2        ;# of coefficients to solve
  53.    x = findgen(nx) # replicate(1., ny)   ;X values at each point
  54.    y = replicate(1.,nx) # findgen(ny)
  55.  
  56.    ut = dblarr(n2, m, /nozero)
  57.    for i=0, degree do for j=0,degree do $
  58.     ut[i*(degree+1) + j, 0] = reform(x^i * y^j, 1, m)
  59.  
  60.    kk = invert(ut # transpose(ut)) # ut
  61.    kx = fltarr(degree+1, degree+1) + float(kk # reform(z, m, 1))
  62.    fit = reform(reform(kx,n2) # ut, nx, ny)
  63.    return, fit
  64. end
  65.  
  66.  
  67.